Identifying complex periodic windows in continuous-time dynamical systems using 

recurrence-based methods 
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The identification of complex periodic windows in the two-dimensional parameter space of certain dynamical 
systems has recently attracted considerable interest. While for discrete systems, a discrimination between pe- 
riodic and chaotic windows can be easily made based on the maximum Lyapunov exponent of the system, this 
remains a challenging task for continuous systems, especially if only short time series are available (e.g., in case 
of experimental data). In this work, we demonstrate that nonlinear measures based on recurrence plots obtained 
from such trajectories provide a practicable alternative for numerically detecting shrimps. Traditional diagonal 
line-based measures of recurrence quantification analysis (RQA) as well as measures from complex network 
theory are shown to allow an excellent classification of periodic and chaotic behavior in parameter space. Using 
the well-studied Rossler system as a benchmark example, we find that the average path length and the clustering 
coefficient of the resulting recurrence networks (RNs) are particularly powerful discriminatory statistics for the 
identification of complex periodic windows. 

PACS numbers: 05.45.Tp, 89.75.Hc, 05.45.Ac 



The investigation of the qualitative behavior in the full 
parameter space of a complex system is a very important, 
but often challenging task. Detailed knowledge about the 
different possible types of dynamics helps understanding 
under which conditions particular states of a system lose 
stability or undergo significant qualitative changes. In 
particular, in experimental settings, the availability of in- 
formation about the underlying patterns in phase space 
allows tuning the critical parameters in such a way that 
one may obtain the desired working conditions. Mathe- 
matically, the corresponding problem is traditionally in- 
vestigated by means of bifurcation theory, which allows 
studying the properties of dynamical transitions in some 
detail QLH]. However, the applicability of available meth- 
ods for identifying bifurcation scenarios and determining 
the parameters at which they take place does often depend 
on the considered system itself. This is especially the case 
when dealing with larger sets of parameters, i.e., operating 
in a two- or even higher-dimensional parameter space, in 
particular for the case of experimental data. In this work, 
we propose some methods based on recurrence properties 
in phase space that allow quantifying dynamically relevant 
properties from available time series, which we harness to 
disentangle the labyrinthine parameter space with respect 
to qualitatively and quantitatively different dynamics. 



I. INTRODUCTION 

The parameter space of nonlinear dynamical systems often 
exhibits a rich variety of qualitatively different types of be- 
havior, such as periodic and chaotic regimes. Especially if 
more than one parameter is responsible for the resulting com- 
plex bifurcation scenario, settings leading to the same types 
of dynamics are often mutually entangled in a rather compli- 



cated way. A specific example are so-called shrimps, a spe- 
cific kind of periodic windows embedded in chaotic regimes 
in the two-dimensional parameter space of a large class of sys- 
tems, which are characterized by a distinct structure consist- 
ing of a main body and four thin legs 0101. The exact proper- 
ties of such structures however depend on the specific system 
under study. In general, shrimps often display self-similarity 
and are regularly organized along some distinguished direc- 
tions . The particular orientation depends on the respective 
stability conditions. When crossing the borders of shrimps in 
different directions, the system typically shows different bi- 
furcation scenarios from periodic dynamics to chaos, e.g., via 
period doubling or via intermittency |4j,|6|] . A detailed stability 
analysis of the periodic solutions contained within shrimps, 
therefore, presents a promising approach for understanding 
the dynamical backbone of these structures. Consequently, 
the study of this particular type of structure has recently at- 
tracted considerable interest. 

The emergence of shrimps has firstly been described in 
great detail for chaotic maps 0B0-@], although correspond- 
ing structures have already been reported in earlier studies 
for both maps and time-continuous systems lflol - [l3ll . In the 
last years, additional efforts have been spent on numerically 
identifying such structures in systems of ordinary differential 
equations (ODEs) as well lllill . Examples include the Rossler 
system 11121 [15141711 . integrate-and-fire models of neurophysi- 
ological oscillations lfl8ll . two coupled parametri cally excited 
oscillators {6], models of laser dynamics [l6l [l9l42lll . the 
damped-driven Duffin g sy stem 12211 . different nonlinear elec- 
tronic circuits lfl^l23fl26tl . chemical oscillators iflrj ITtI 27], 
the Lorenz-84 low-order atmospheric circulation model 11611 . 
or a vibro-impact system 11281. Recently, shrimp structures 
have been experimentally found in the chaotic Chua ||29[ l30ll 
and Nishio-Inaba circuits l3ll) . which underlines the increas- 
ing importance of proper numerical algorithms for automat- 
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ically identifying such periodic windows in both theoretical 
and experimental studies, possibly even in the presence of 
noise. 

Typically, at the boundaries of a shrimp, small inaccuracies 
of the parameters induce dramatic changes in the resulting dy- 
namics. Therefore, these structures can hardly be uncovered 
by existing analytical methods based on linear stability analy- 
sis, impossible in particular at the boundaries. In most recent 
works for continuous systems, the maximum Lyapunov ex- 
ponent has been estimated numerically, which allows dis- 
tinguishing periodic and chaotic dynamics since Ai = for 
periodic orbits and Ai > for chaotic ones. Recently, an al- 
ternative method has been proposed for uncovering shrimps 
in systems described by ODEs 03] based on the concept of 
correlation entropy A'2, which can be considered as a lower 
bound for the sum of all positive Lyapunov exponents of the 
system ll32ll . It has been demonstrated that considering K2 in- 
deed yields results of comparable accuracy as the maximum 
Lyapunov exponent method [6] and is more efficient, because 
less data points are needed. 

A convenient way for numerically estimating A'2 is using 
recurrence plots (RPs), an efficient technique of nonlinear 
time series analysis. Given a trajectory of a dynamical system 
consisting of different values x,, where i indicates the time of 
observation, the corresponding RP is defined as 1I33I l34tl 

JZ iJ (e) = e(e-||x i -xj||) ) (1) 

where O(-) is the Heaviside function, ||xj — Xj|[ = dij is the 
distance between two observations and Xj in phase space 
(which will be measured in terms of the maximum norm in 
the following), and e a pre-defined threshold for the prox- 
imity of two states in phase space, i.e., for distinguishing 
whether or not two observations are neighbors in phase space. 
Hence, the basic mathematical structure describing a RP is 
the binary recurrence matrix Ri.j- Visualizing this matrix by 
black (Rij = 1) and white {Rij = 0) dots, different types 
of dynamics can be identified in terms of different features 
of line structures (including discrete points, blocks, and ex- 
tended diagonal or vertical lines), which can be quantitatively 
assessed in terms of recurrence quantification analysis (RQA, 
see Sec. Ill At . 

The recurrence plot based estimation of A'2 involves two 
main steps: (i) computing the cumulative probability distribu- 
tions of the lengths of diagonal lines, and (ii) properly select- 
ing a scaling region in dependence on the diagonal line length 
I and evaluating its characteristic parameters |35 1. Since the 
second step assumes the existence of sufficient convergence in 
a reasonable regime, in practical applications various values of 
e need to be considered for properly detecting the correspond- 
ing scaling region. Hence, this approach is partially subjective 
and depends on the specific choice of the scaling region. 

In this work, we suggest and thoroughly study alternative 
approaches to analyze quantitatively the corresponding prop- 
erties of the recurrence matrix. More specifically, we ap- 
ply measures from complex network theory to the recurrence 
structures and compare their performance with that of more 
traditional RQA measures. The corresponding framework of 
recurrence networks (RNs), i.e., the idea of interpreting the 



recurrence matrix Ri j of Eq. (Q]i as the adjacency matrix Aij 
of an undirected, unweighted complex network by setting 

where dij is the Kronecker delta, has been recently proposed 
for identifying dynamical transitions in model systems such 
as the logistic map l36ll . as well as detecting subtle changes 
in a marine dust flow record. It should be noted that simi- 
lar approaches of understanding the properties of time series 
from complex network perspectives have been suggested in 
parallel by various authors H [H (see JH |H for further 
references). 

While most classical RQA measures are based on line struc- 
tures in the RP, complex network measures reveal higher- 
order properties of the phase space density of states of the con- 
sidered dynamical system 113911 (see Sec. Ill Bb . However, the 
question whether the complex network approach is more suit- 
able than conventional RQA for distinguishing different types 
of dynamics has not yet been systematically addressed. In 
this work, we provide a corresponding comparative analysis 
using some particularly well-suited measures of both types. 
Specifically, we demonstrate clear advantages of complex net- 
work measures in identifying shrimp structures or, even more, 
arbitrary complex periodic windows in the two-dimensional 
parameter space of certain dynamical systems by taking the 
well-studied Rossler system as an illustrative benchmark ex- 
ample. The performance of different measures is compared 
by means of several statistical tests based on the resulting pat- 
terns in the parameter space. 

The remainder of this paper is organized as follows: In 
Sec. |nl we briefly review some RQA as well as complex net- 
work measures and outline important technical issues arising 
when using these measures, in particular, the dependence on 
the choice of the recurrence threshold e and other parame- 
ters. A specific part of the two-dimensional parameter space 
of the Rossler system that is known to contain shrimps is ex- 
amined in Sec.|nl]by means of the maximum Lyapunov expo- 
nent Ai, which serves as a reference for our results obtained 
using recurrence-based methods in Sec. [IV] The main conclu- 
sions of our work are summarized in Sec.[Vl 

II. METHODS 
A. Recurrence quantification analysis (RQA) 

Recurrence quantification analysis (RQA) has been intro- 
duced fo r q uantifying the presence of specific structures in 
RPs Il4ll444ll . During the last years, numerous methodologi- 
cal developments and applications in various fields of science 
have been reported (for a corresponding review, see 1341 ). 
Most traditional RQA measures are based on the length dis- 
tributions of diagonal or vertical lines. In this work, we will 
particularly make use of the following three quantities: 

• The recurrence rate 

2 x 

RR= N ( N ^i)2^ R ^> (3) 
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measures the fraction of recurrence points in a RP and, 
hence, gives the mean probability of recurrences in the 
system. 

• The determinism DET is defined as the percentage of 
recurrence points belonging to diagonal lines of at least 
length l min (see Sec. IIII Cl for details). 



DET 



N 

I=i„ 



W(l) 



(4) 



where P(l) denotes the probability of finding a diagonal 
line of length I in the RP. This measure quantifies the 
predictability of a system. 

• The average diagonal line length Lmean, defined as 



L 



MEAN 



(5) 



characterizes the average time that two segments of a 
trajectory stay in the vicinity of each other, and is re- 
lated to the mean predictability time. 

In addition to (diagonal as well as vertical) line-based RQA 
measures, in some cases, parameters characterizing recur- 
rence times (i.e., vertical distances between two recurrence 
points) based on the RPs have been suggested to be suitable 
for quantifying basic properties of the recorded dynamics l45l - 

Let us now consider in what sense RQA measures behave 
differently in the presence of periodic and chaotic dynamics. 
In a RP, a periodic orbit is reflected by long non-interrupted 
diagonal lines that are separated by a constant offset, which 
corresponds to the period of the oscillation (Fig.[T]\). In con- 
trast, for a chaotic trajectory, the distances between diagonal 
lines are not constant due to multiple time scales present in the 
system (Fig.QJ}). Furthermore, the existing diagonals are in- 
terrupted because of the exponential divergence of nearby tra- 
jectories. Therefore, we expect that both DET and Lmean 
typically have larger values for a periodic trajectory than for a 
chaotic one. 
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FIG. 1: Recurrence plots for a periodic (A) and a chaotic (B) trajec- 
tory of the Rossler system ((9} (see Sec. lIIIBl for details). 



B. Quantitative analysis of recurrence networks 

Recurrence networks (RNs) are spatial networks represent- 
ing local neighborhood relationships in the phase space of 
a dynamical system. In a RN, every sampled point is rep- 
resented by a vertex v, whereas edges indicate pairs of ob- 
servations whose mutual distance in phase space is smaller 
than the pre-defined threshold e. The resulting edge density 
p, i.e., the relative fraction of edges that are actually present 
in comparison to a fully connected network, is then given by 
RR (see Eq. (O). Hence, characterizing topological prop- 
erties of RPs by sophisticated network-theoretic quantifiers 
allows retrieving additional information about higher-order 
phase space properties of the system l39ll . It has been shown 
that these measures are able to capture dynamical transitions 
in complex systems, such as in the logistic map with changing 
control parameter or a real-world paleoclimatic time series, 
demonstrating that network-theoretic features provide addi- 
tional and complementary information when compared to tra- 
ditional RQA measures Pol . 

In this work, we particularly consider the global cluster- 
ing coefficient C and the average path length £ of a RN. The 
application of other measures (e.g., betweenness centrality) 
is straightforward, but we argue that the description of net- 
work topology by a scalar parameter instead of a distribution 
of vertex- or edge-based statistics might be more suitable for 
detecting and quantifying qualitative changes in the dynamics 
of a system. Moreover, one has to note that state-of-the-art 
numerical algorithms for computing other complex network 
measures often require considerably larger computational ef- 
forts than estimating traditional RQA measures. 

The local clustering coefficient C v characterizes the phase 
space geometry of an attractor in the e-neighborhood of a ver- 
tex v. Specifically, C v gives the probability that two randomly 
chosen neighbors of v are also neighbors Bill , i.e., 



r 



k v {k v 1) 



N: 



(6) 



where k v is the degree centrality (i.e., the number of neighbors 
of v, which coincides with the local recurrence rate RR V ), and 
is the total number of closed triangular subgraphs includ- 
ing v, which is normalized by the maximum possible value 
k v {k v — l)/2. For vertices of degree k v = or 1 (i.e., isolated 
or tree-like vertices, respectively), C v = by definition since 
such vertices cannot participate in triangles. Instead of study- 
ing C v individually for all vertices of a RN, we consider its 
average value taken over all vertices of a network, the global 
clustering coefficient 



1 N 



(7) 



as a global characteristic parameter of network topology. 

The local clustering coefficient C v is related to the ef- 
fective dimensionality of the set of observations in the e- 
neighborhood of a vertex v l39ll . In particular, we find specific 
dependencies for both continuous and discrete dynamical sys- 
tems: 
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(i) For discrete systems, periodic orbits consist only of a fi- 
nite set of points. This implies that for e <C A with 
A > maxij di.j being the attractor diameter, the RN is 
decomposed into disjoint components with multiple ver- 
tices being located at the same point in phase space. As 
a consequence, the individual components are fully con- 
nected for all e > 0, which leads to C = 1. In contrast, 
for chaotic trajectories, sufficiently small e« A results 
inC < 1 & 

(ii) For a continuous system, it is a well-established fact that 
if a chaotic trajectory enters the neighborhood of an un- 
stable periodic orbit (UPO), it stays within this neighbor- 
hood for a certain time 14911 . As a consequence, states 
accumulate along this UPO instead of homogeneously 
filling the phase space in the corresponding neighbor- 
hood (in particular if we consider UPOs of lower period). 
Since from the theory of spatial random graphs 115011 
it is known that the average clustering coefficient in- 
creases with decreasing spatial dimension of the space 
in which a network is embedded, it follows that parts of 
a chaotic attractor that are close to (low-periodic) UPOs 
are characterized by higher values of C v in the resulting 
RNs ll39ll . Following a related argument, we expect that 
states on a periodic orbit of a continuous system are also 
characterized by high values of C v , which implies high 
C. In general, however, for chaotic trajectories the filling 
of the phase space with observed states is more homoge- 
neous than for periodic ones, which leads to a tendency 
towards lower values of C v and, hence, C. 

A second global measure of network topology is the av- 
erage path length. Since we consider RNs as undirected 
and unweighted, we define all edges to be of unit length in 
terms of graph (geodesic) distance. The graph distance be- 
tween any two vertices of the network is defined as the length 
of the shortest path between them. Hence, the shortest path 
length li.j in the RN gives the minimum number of edges that 
have to be passed on a graph between two vertices i and j. 
Accordingly, kj is related to the phase space distance of the 
associated states, but not to the temporal evolution of the sys- 
tem (Fig. 

The average path length C is defined as the mean value of 
the shortest path lengths lij taken over all pairs of vertices 

For disconnected pairs of vertices, the shortest path length is 
set to zero by definition. Note that in most practical applica- 
tions, this has no major impact on the corresponding statistics. 
The average path length is related to the average separation of 
states in phase space, which measures the size of the attrac- 
tor in units of e. Note that since metric distances in the phase 
space are conserved, this imp lies an approximately recipro- 
cal relationship C oc [39], with the exception of periodic 
orbits of discrete maps (see the arguments given above). 



For the behavior of the average path length C, one has to 
carefully distinguish between continuous and discrete systems 
again: 

(i) For discrete systems, the structure of a periodic orbit in 
phase space implies C = 1 by definition for e <C A l36ll . 
In turn, for chaotic trajectories, there is a continuum of 
possible values, so that C > 1 for e <C A. Hence, the 
average path length of a periodic orbit is smaller than 
that of a chaotic one. 

(ii) For continuous systems, however, a periodic trajectory 
has a completely different structure in phase space, 
which means that C is approximately determined by the 
curve length of the orbit in the phase space (Fig.|2j^)- In 
contrast, for chaotic trajectories, there are "short-cuts" 
that allow reaching one position in phase space from an- 
other with a lower number of steps than by following the 
trajectory l39ll (Fig.|2j3). Consequently, for a comparable 
attractor diameter A and the same recurrence threshold 
e, C can be expected to have lower values for chaotic or- 
bits. Alternatively, it is desirable to fix RR instead of e 
to obtain RNs with approximately the same numbers of 
edges. A detailed discussion on the advantages of this 
approach is presented in Sec. MIDI In this case, we note 
that the different geometric structure of the attractor in 
both the periodic and chaotic regime even enhances the 
effect of shortcuts, so that under general conditions, pe- 
riodic orbits of continuous systems are characterized by 
larger £ than chaotic ones. 




-80 -40 40 80 -80 -40 40 80 



X X 

FIG. 2: (Color online) Illustration of the shortest path length for two 
example trajectories of the Rossler system l[9} in (A) periodic and 
(B) chaotic regime. The square is a schematic projection of the re- 
currence neighborhood to the (x,y) plane. In these two particular 
examples, if^j = 24 and = 19 for e = 7.0 (maximum norm). 

Concerning the behavior of the average path length as 
shown in Fig. [2] we have to point out that differences between 
periodic and chaotic dynamics require that transient behav- 
ior (e.g., during some initial phase of the dynamics before the 
attractor has been reached by the considered trajectory) has 
been sufficiently removed, since such transients could also 
result in artificial short-cuts in phase space for periodic sys- 
tems. While this problem can be easily solved in numerical 
studies by simply discarding the transients, it can contribute 
additional uncertainties in experimental studies where only a 
limited amount of data is available. In a similar way, using 
C for discriminating between periodic dynamics and chaos in 
noisy experimental time series is expected to lead to the same 
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problem. However, in such cases, other network-theoretic 
quantities (e.g., C) still provide feasible alternatives, to which 
the mentioned conceptual problem does not apply (at least for 
sufficiently small noise levels that do not exceed the order of 
magnitude of the recurrence threshold e). In addition, there 
are other recurrence-based techniques such as recurrence time 
statistics, which allow distinguishing regular and chaotic dy- 
namics even in challenging situations, e.g., quasi-periodicity 
versus weakly chaotic behavior ll45l Bill . 

In summary, C shows a consistent behavior for discrete and 
continuous systems, whereas C performs inversely with re- 
spect to periodic and chaotic dynamics. This implies, that C 
cannot be used alone to classify the dynamics of time series 
from real systems, where the nature of the underlying pro- 
cess is unknown. Our results presented in this work and else- 
where j3(J suggest, however, that C is nevertheless very well 
able to distinguish between periodic and chaotic behavior in 
the parameter space of the same dynamical system. 

HI. MODEL: CHAOTIC ROSSLER SYSTEM 

A. Two-dimensional parameter space 

As an illustrative example for a time-continuous dynamical 
system that shows shrimp structures in its parameter space, we 
study the Rossler system 

/ dx dy dz\ , . .. 

\~dt''~dt , 'dt) = ~ Z,X + ay + X ~ C '' ' 

In this system, c is often regarded as the standard bifurcation 
parameter, whereas a and b mainly change the shape of the 
attractor. In this work, however, we vary a = b on the one 
hand, and c on the other hand, which yields a two-dimensional 
parameter space of (c, a). One of the parameter regions of 
special interest is (c, a) £ [20,45] x [0.2,0.3], where the 
chaotic regime is intermingled with periodic windows in a 
complex manner lfl5ll . Thus, in the following, we pay spe- 
cial attention to this part of the parameter space. Note that 
shrimp structures are not confined to this particular plane in 
the three-dimensional parameter space of the Rossler system 
and, hence, investigating a differently oriented plane as in lfl6ll 
would not qualitatively alter the results of our study. 

For further analysis, the above mentioned part of the pa- 
rameter space of (c, a) is divided into 1,000 x 1,000 grid 
points, which results in the step size 0.0001 in a £ [0.2, 0.3] 
and 0.025 in c 6 [20,45]. We first consider the maximum 
Lyapunov exponent of the resulting systems as a character- 
istic measure to determine whether the trajectories resulting 
from each parameter combination are chaotic or periodic. For 
this purpose, all three Lyapunov exponents A, have been com- 
puted based on Eq. (0 by numerically solving the associated 
variational equations |52J. Numerical integration is carried 
out using a fourth-order Runge-Kutta integrator with a fixed 
step size of h = 0.001 time units and randomly chosen initial 
conditions (xo,yo,zo) S [0,1] x [0,1] x [0,1]. In order to 
provide sufficient data for an accurate estimation of Lyapunov 



exponents, the integration is performed for sufficiently long 
time (until t = 10, 000, i.e., involving N = 10, 000, 000 data 
points with a sampling At corresponding to the considered 
integration step h). For computing the different recurrence- 
based measures in a second step of analysis, the initial tran- 
sients (cf. Sec. Ill Bt have been discarded by removing the first 
500, 000 iterations (i.e., all data until t = 500) from the simu- 
lated trajectories. Moreover, for all further calculations based 
on the recurrence matrices, a coarser sampling of At = 0.2 
(i.e., 200 integration steps) has been considered. For each tra- 
jectory, the sampling has been terminated after N = 5, 000 
points (corresponding to about 150 ~ 250 oscillations of the 
system) have been obtained, yielding much shorter time se- 
ries than those used for estimating Lyapunov exponents. Note 
that the specific choice of sampling time has a certain influ- 
ence on the results discussed in the following. Specifically, 
the sampling corresponding to the optimum resolution of dy- 
namically relevant features is expected to vary among differ- 
ent parameter settings. In this respect, the considered value of 
At represents a reasonable and practically tractable choice. 




FIG. 3: (Color online) Maximum Lyapunov exponent Ai in the (c, a) 
parameter plane of the Rossler system (|9). Regions with Ai = indi- 
cate periodic dynamics, those with large Ai correspond to a strongly 
chaotic behavior. Asterisks indicate the parameter combinations used 
as examples in Sec.Hllandlllll 



The resulting maximum Lyapunov exponent Ai reveals a 
rich structure of chaotic and periodic dynamics in the param- 
eter space (Fig. 0. Inside the chaotic regions, several well 
pronounced periodic bands are identified. Moreover, in the 
center of the plot, a special periodic window of interest is 
found, with structures resembling a head and four main thin 
legs. These particular swallow-like structures can be identi- 
fied as shrimps lfl5ll . Additional zooms into the parameter 
space uncover numerous fine structures, which correspond to 
secondary shrimps. 
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B. Prototypical trajectories 

We first illustrate the potentials of recurrence-based ap- 
proaches for discriminating between periodic and chaotic 
dynamics by studying trajectories obtained from two dis- 
tinct parameter combinations. Based on the parame- 
ter space shown in Fig. [3] the dynamics is periodic 
for {a,b,c) = (0.245,0.245,35), while chaotic for 
(a,b,c) = (0.29,0.29,35). The Lyapunov exponents 
are (0, —0.223, —32.683) for the periodic trajectory, and 
(0.151,0, —32.517) for the chaotic one. The typical phase 
space distances of states on the two trajectories are of the same 
order of magnitude (Tab. Q}, which is illustrated by the prob- 
ability distribution of the mutual distances di.j in Fig. [4j\. 
However, the detailed relationship of RR(e) is clearly dif- 
ferent: for a low (but fixed) RR, we have to consider much 
higher e for chaotic trajectories than for periodic ones, while 
for larger RR(e), the opposite behavior is found (Fig. |4j3). 
RR corresponds to the correlation sum, the slope of which 
in a double logarithmic plot gives the correlation dimension 
£>2- Therefore, RR(e) having a larger slope for the chaotic 
trajectory than for the periodic one (Fig.|4j3) indicates that the 
correlation dimension of the former is higher than that of the 
latter. 





periodic 


chaotic 


Ai 


0.0 ± 0.0003 


0.15 ±0.0002 


(dij) 


58.95 ±0.14 


57.42 ± 3.03 




476.69 ±0.15 


575.12 ±23.99 


e (RR = 0.02) 


4.79 ±0.14 


6.50 ±0.53 


DET 


0.99 ± 0.01 


0.97 ±0.01 


Lmean 


18.51 ±0.18 


10.60 ± 0.34 


C 


0.77 ±0.001 


0.62 ±0.004 


c 


45.30 ± 1.83 


9.19 ±0.37 



TABLE I: Maximum Lyapunov exponents Ai (N = 10, 000, 000, 
At — 0.001), mean and maximum separation of points in phase 
space (N = 5, 000, At = 0.2) and resulting recurrence threshold 
e (maximum norm) for RR = 0.02, and RQA (l m in = 2) and net- 
work measures for two parameter combinations (see text), represent- 
ing periodic and chaotic regimes of the Rossler system. The error 
bars correspond to the standard deviation obtained from 100 realiza- 
tions with different initial conditions. Note that the large variance of 
the metric quantities e and dij for the chaotic trajectory is a com- 
mon result when working with short time series and different initial 
conditions 15311 . 

Both RQA and complex network measures highlight dif- 
ferences in the topological structure between the periodic and 
chaotic regimes (Tab. JJ. We observe that for the considered 
recurrence rate (RR = 0.02), DET, Lmean, C and C have 
lower values for chaotic trajectories than for periodic ones 
(which distinctively differs from the behavior of the maximum 
Lyapunov exponent Ai), while the difference is significant for 
the latter three measures. However, before being able to gen- 
eralize these findings, the robustness of the aforementioned 
features has to be critically examined. Our corresponding re- 
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FIG. 4: (A) Probability distributions p(dij) of mutual distances di,j 
(maximum norm) between states on one realization (N = 5, 000) of 
periodic (solid) and chaotic (dashed) trajectories, respectively (see 
text). (B) Dependence of the recurrence rate RR on the recurrence 
threshold e for a periodic and chaotic regime. The error bars corre- 
spond to the standard deviation obtained from 100 realizations with 
different initial conditions. 



suits are summarized in the following two subsections. 



C. Dependence on l min 

For the computation of most line-based RQA measures 
(with some exceptions, such as RR or the maximum diag- 
onal line length), one has to choose a proper minimal line 
length lmin, to avoid a bias due to oversampling l34ll . The 
particular values of these measures do significantly depend on 
lmin (Fig. |5]), although the discrepancy between the values of 
DET and Lmean obtained for periodic and chaotic trajecto- 
ries remains qualitatively unchanged (in general, the values of 
DET and Lm ean are larger for a periodic trajectory than for 
a chaotic one). Note that DET and Lmean are more robust 
against noise than other measures like the maximum diagonal 
line length jH]. 
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FIG. 5: Dependence of the RQA measures DET (A) and Lmean 
(B) on l m in for periodic (solid) and chaotic (dashed) trajectories. 
The error bars indicate the standard deviation obtained from 100 re- 
alizations of the Rossler system ® with TV = 5, 000, At = 0.2, 
RR — 0.02, and different initial conditions. 



D. Dependence on the recurrence threshold 

Since all considered measures are based on the recurrence 
matrix, their values necessarily depend on the choice of the 
recurrence threshold e. So far, there is no universal threshold 
selection criterion for the RP computation. On the one hand, 
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if e is chosen too small, there are almost no recurrence points 
and, hence, no feasible information on the recurrence structure 
of the system. On the other hand, if e is too large, almost 
every point is a neighbor of every other point, which leads 
to numerous artifacts. Hence, we have to seek a compromise 
for choosing a reasonable value of e. One rule of thumb is 
choosing e in such a way that it lies in the scaling regime of 
the double logarithmic plot of the correlation integral versus 
e. This rule coincides with the classical strategy for estimating 
the correlation dimension D 2 using the Grassberger-Procaccia 
algorithm ifsill . Following independent arguments, Schinkel 
et al. f55ll suggested choosing e corresponding to at most 5% 
of the maximum attractor diameter in phase space. In a similar 
way, RR < 0.05 has been found a reasonable choice in the 
analysis of RNs SIM El- 

The dependence of RQA and RN measures on e is shown 
in Fig. [6] One observes that all previously discussed measures 
are clearly able to distinguish between periodic and chaotic 
dynamics, i.e., they show significant differences in a broad 
interval of e (corresponding to mean recurrence rates of about 
l%to5%). 



A 0.99r 



A 0.99 



0.96 



C 0.8 



U0.7 



H-i-H-H+H-H-I 



B 20r 



12 



M-l-H-H-I-H-H-I 

456789 10 456789 10 

£ 8 



4 5 6 7 8 9 10 




0.98 



0.97 f- 



0.96 L 



B 19 



0.01 0.02 0.03 0.04 0.05 
RR 




0.01 0.02 0.03 0.04 0.05 
RR 



C 0.8 r 
0.75 

iO 0.7 
0.65 
0.6 



0.01 0.02 0.03 0.04 0.05 
RR 




0.01 0.02 0.03 0.04 0.05 
RR 



FIG. 7: Same as in Fig. |6]for the dependence of these measures on 
the recurrence rate RR. 



the numerical algorithms for estimating all RP-based mea- 
sures, fixing e would lead to values of RR and other RQA 
measures that cannot be compared in a meaningful way. A re- 
lated approach using the dependence on RR has already been 
applied for an automatized detection of the scaling region in 
the recurrence-based estimation of K2 for the study of a large 
scale system lf5r3l . 

According to the above arguments, we will use a fixed value 
of RR in all following calculations. Because of this, in a first 
step of analysis, we study the dependence of the different mea- 
sures on RR in some detail similar to the dependence on e 
discussed above. Figure|7]reveals that the difference between 
periodic and chaotic orbits remains qualitatively unchanged 
and does not depend on the specific choice of RR. Specifi- 
cally, in all four cases, periodic trajectories are characterized 
by larger values of the different measures. 
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FIG. 6: Dependence of the RQA measures DET (A), Lmean (B), 
and the network measures C (C) and C (D) on the recurrence thresh- 
old e for periodic (solid) and chaotic (dashed) trajectories (see text). 
The error bars indicate the standard deviation obtained from 100 re- 
alizations of the Rossler system {9} with N = 5, 000, At = 0.2, 
Imin = 2, and different initial conditions. 

Instead of fixing the recurrence threshold e, it may be desir- 
able to compare different situations by using RPs with a fixed 
value of RR. Firstly, the resulting RNs have approximately 
the same number of edges, which allows comparing the re- 
sulting topological properties of different networks more ob- 
jectively. Secondly, it has been shown in Fig. 2(3 that there is 
a crossover between the RR of periodic and chaotic trajecto- 
ries obtained with the same e, which is related to the fact that 
RR(e) (which corresponds to the correlation sum) increases 
exponentially for a chaotic trajectory llsill . Finally, we have 
to note that the amplitudes of the (chaotic or periodic) oscilla- 
tions of the Rossler system change with the specifically cho- 
sen control parameters of the system. In order to automatize 



IV. RESULTS 

In the following, we study the performance of the different 
recurrence-based measures for identifying shrimp structures 
in the Rossler system. For this purpose, we take Ai as a refer- 
ence, since this measure per definition discriminates between 
periodic (Ai = 0) and chaotic (Ai > 0) dynamics. Note that 
in order to reach reliable estimates of Ai, typically a large 
amount of data is required. In contrast, recent results on RQA 
and RN measures 13611 suggest that these methods allow iden- 
tifying dynamical transitions in complex systems using much 
shorter time series. 

In order to systematically test the applicability of 
recurrence-based methods for identifying shrimps, we use the 
same grid with 1, 000 x 1, 000 pairs of points in the (c, a)- 
parameter plane as in Fig. [3] For each set of parameters, we 
consider time series of A = 5, 000 data points sampled with 
a fixed time step of At = 0.2 for estimating all RQA and RN 
measures. Initial transients have been removed from the data 
as described in Sec.lIIIK. We stress that measures originated 
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from both RQA and network theory are computed from the 
same recurrence matrices R^j (which have been constructed 
from the original three-dimensional coordinates of the system 
without embedding), so that the resulting structures in the pa- 
rameter space are well comparable. In the following, we will 
consistently use l m in = 2 and RR = 0.02. 

A. Behavior of individual measures 

The values of the four chosen measures in dependence on 
the parameters a and c are shown in Fig. [8] We find that all 
four measures are indeed able to identify periodic windows, in 
particular, those associated with shrimp structures. However, 
the discriminatory power of the individual parameters for dis- 
tinguishing between periodic and chaotic regions is notably 
different. Specifically, for the chosen values of l m i n and RR 
the contrast in the values of DET obtained for both regimes 
is relatively weak (see Figs. [5]and[7]and Tab. JJ, whereas the 
other three measures show a much larger range of values. 

Unlike DET and Lmean, C and C resolve all peri- 
odic regimes including band structures and even secondary 
shrimps. The main remaining difference to the structures ob- 
tained by using the maximum Lyapunov exponent (Fig.O is 
found in the broad periodic band in the upper right corner of 
the parameter plane. In this region of parameter space, the 
system shows a complex bifurcation scenario, including dif- 
ferent routes from periodic behavior to chaos similar to those 
associated with shrimp structures in other continuous dynam- 
ical systems IQj. These differences are partially related to nu- 
merical inaccuracies in the presence of time series with a finite 
length and become less pronounced if longer simulations are 
used. Since no analytical solutions of the system are avail- 
able for this specific region, additional application of comple- 
mentary numerical methods (e.g., Lyapunov exponents, return 
maps, Poincare sections, etc.) would be necessary for com- 
pletely identifying the associated bifurcation scenario, which 
is however out of the scope of this paper. 

B. Correlations between recurrence-based measures and Ai 

We now perform a more detailed quantitative analysis of 
the power of the individual measures as discriminatory statis- 
tics between periodic and chaotic dynamics. For this pur- 
pose, the differences between the structures obtained using 
the discussed measures and those revealed by Ai can be in- 
terpreted as a lack of performance. Since the different mea- 
sures are characterized by strongly differing probability distri- 
bution functions (PDFs), we seek a statistics that allows prop- 
erly quantifying deviations between the respective structures 
in the parameter space. The first insight into this question is 
provided by the correlation coefficients between Ai and the 
other measures, which are in all cases clearly significant, but 
with a negative sign (Tab. Hill. This result supports the findings 
in Sec.|nI]for the two example trajectories. 

In order to study how strongly the patterns of the differ- 
ent measures in the (c, a)-plane resemble those found for the 



maximum Lyapunov exponent, we investigate the point-wise 
difference of the corresponding cumulative distribution func- 
tions (CDFs), i.e., 

AP(X 1 ,x) = P(X 1 )-P{x), (10) 

where P{x) is (for a given (c, a)-combination) the empirical 
value of the CDF of the measure x obtained from all 1 ,000,000 
parameter combinations. In order to simplify the notation, the 
arguments (c, a) of AP(Xi,x) will be suppressed. Note that 
the maximum absolute value of AP, commonly denoted D, 
corresponds to the Kolmogorov-Smirnov statistics frequently 
used for testing the homogeneity of the probability distribu- 
tions of two samples ll57tl . Since all four recurrence -based 
measures are negatively correlated with Ai, we will identify x 
with 1 — DET, 1 — Lmean, 1 — C, and 1 — L, respectively. 





P 


„2 


DET 


-0.75 


0.21 


Lmean 


-0.81 


0.18 


C 


-0.70 


0.23 


C 


-0.66 


0.24 



TABLE II: Overall performance indicators obtained from a point- 
wise comparison of the values of the maximum Lyapunov exponent 
Ai and the different RQA and network measures: Spearman's p and 
the standard deviation a 2 of the CDF differences AP(Ai, x). For 
simplicity, the arguments of the different characteristics have been 
omitted. 

Due to the significant correlations between the differ- 
ent measures (Tab. UD, we expect that the CDF difference 
AP(Ai,a;) is close to zero in large parts of the parameter 
space. The patterns of the CDF differences are shown in 
Fig. [9] In all cases, values close to zero can often be observed 
in periodic windows (although differences appear especially 
in the secondary shrimps), whereas the shortness of the con- 
sidered time series leads to significant deviations from zero in 
the chaotic regions. The standard deviations a\ P of the CDF 
field provide a rough impression of the amount of incorrectly 
identified regimes (partially due to the finite length of time 
series), which will be quantitatively characterized in the next 
section. 



C. Probability of classification errors 

The general treatment in Sec. |IV B| do not yet allow so- 
phisticated conclusions on which particular measure is most 
appropriate for detecting shrimps in continuous systems. In 
order to assess the discriminatory power of all four measures 
in more detail, we subject the resulting patterns in parameter 
space (Fig. [8J to further statistical analysis. For this purpose, 
we explicitly make use of the fact that the (c, a) -parameter 
plane of the Rossler system is composed of two sets of param- 
eter combinations belonging to trajectories with periodic and 
chaotic dynamics, respectively. However, the corresponding 
group structure is not exactly known, since there are numerous 
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FIG. 8: (Color online) RQA measures DET (A), Lmean (B), and network measures C (C) and C (D) in the (c, a) parameter plane of the 
Rossler system (|9). 



(c, a) -combinations which lead to small values of Ai (Fig. [Tot. 
These corresponding values could hence be related to either 
weakly chaotic behavior or periodic dynamics that cannot be 
exactly detected due to the numerical precision. Note that no 
fixed points (Ai < 0) are found in our parameter space. In or- 
der to obtain a robust approximation of the two distinct groups 
of parameter combinations leading to Ai = and Ai > 0, re- 
spectively, we refer to a (variable) critical value A* (Fig. ITOt 
for defining two disjoint sets 



Si(X") := {(c,o)|Ai(c,a) < A*} 
S 2 (A*) := {(c, a)|Ai(c, a) > A*} 



with group sizes n\ and 77,2 
case). 



m (n = 1,000, 000 in our 



1. F and U tests 

One possibility for assessing the discriminatory power of 
the different recurrence-based measures is taking the respec- 
tive groups Si and S2 (for different values of A* in a reason- 
able range, i.e., A* e [0,0.05], see Fig. ITOb and statistically 
quantifying whether or not main statistical characteristics of 
the distributions p(x\Si) of the different measures x obtained 
for both groups Si and S2 differ significantly. For one spe- 
cific choice of A*, these distributions are shown in Fig. QT| 
Formally, this question corresponds to a one-way analysis of 
variance (ANOVA) problem I158ll . with the factor being deter- 
mined by two classes of values of Ai. In order to evaluate 
whether the group means do significantly differ (in compar- 
ison with the respecti ve g roup variances), the F-test is used 
with the test statistics 115911 



T71 (Ml 

t = mri2 — 2 



M2) 2 = 2 

n — t 1 



msi + T12S2 



(11) 



10 



A 0.3 



0.28 



0.26 



05 



0.24 



0.22 



0.2 



C 0.3 



0.28 



0.26 



03 



0.24 



0.22 



AP(X, 1 , DET) 





\ 



































w 



0.97 B 0.3 1 



0.28 



0.47 



0.26 



-0.02 03 



i 



-0.51 



0.24 



0.22 



20 25 30 35 40 45 

c 



-1.00 0.2 1 



£P(Ai, C) 




0.99 D 0.3 



0.28 



0.49 



0.26 



-0.00 03 



-0.50 



0.24 



0.22 



ap <V l mean) 




\ 



20 25 30 35 40 45 

c 



AP(Ai, C) 



-0.99 0.2 




I 

n 



0.91 



0.44 



-0.04 



-0.52 



-0.99 



0.98 



0.48 



-0.01 



-0.50 



'-1.00 



FIG. 9: (Color online) CDF differences AP(Ai, ■) between the maximum Lyapunov exponent Ai and the RQA measures DET (A), Lmean 
(B), and network measures C (C) and £ (D) in the (c, a) parameter plane of the Rossler system l[9j. 



where t is the test statistics of a corresponding i-test |6C 
Since we are aware that the values of F (or, alternatively, t) 
may be misleading if the underlying sample PDFs are strongly 
non-Gaussian (see Fig. [TTJ, our results are complemented by 
those of a corresponding distribution-free test. Specifically, 
we compute the value of the test statistics U of the Mann- 
Whitney U -test i6ll l62ll against equality of the medians of 
two distributions, which can be considered as the equivalent 
of an J 71 - test performed on the sets of rank numbers. 

For both tests, we find that for almost all choices of A*, C 
and C show the highest values of the respective test statistics, 
indicating that the discriminatory power for distinguishing be- 
tween both sets is the largest for these measures (Fig.[T2K,B). 
Note that in all cases, the probability values for rejecting the 
respective null hypotheses (i.e., equality of group means and 
medians, respectively) are close to 100% due to the large sam- 



ple size. The results obtained using both test statistics are sup- 
ported by a numerical approximation of the associated overlap 
integral 



i(x) 



dxp(x\Si) p(x\S 2 ) 



(12) 



of the PDFs of the recurrence based measures x for both 
groups (Fig. [TTJ, the values of which are shown in Fig.fTZC. 

The advantages of RN measures (at least for this particu- 
lar case) become particularly apparent when visually inspect- 
ing Fig. [TTJ The overlap of the PDFs p(x\Si) can be seen to 
be substantially smaller for the network quantifiers C and C 
(Fig.fTTC.D). than for the RQA measures DET and L M ean 
(Fig.[ll&B). 
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FIG. 10: Probability distribution function of the maximum Lyapunov 
exponent Ai obtained from all 1,000,000 parameter combinations in 
the considered (c, a) plane of the Rossler system ©. 
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FIG. 12: Measures for the discriminatory skills of the different 
recurrence-based measures DET (•), Lmean (P), C (+), and 
C (>), obtained from a comparison with the results derived using 
the maximum Lyapunov exponent in dependence on the choice of 
A*: (A) t-test statistics, (B) t/-test statistics, (C) overlap integral ^ 
(Eq. d!2t). and (D) relative frequency of false detections p/, using 
the same quantiles of Ai and the respective measures. 



frequency of "grouping errors" in comparison to Ai, followed 
by C, DET and L MEAN . 



3. ROC analysis 



FIG. 11: PDFs p(x\Si) of RQA and RN measures for parameter 
combinations (c, a) yielding maximum Lyapunov exponents Ai < 
A* (Si) and Ai > A* (S2) with A* = 0.01, respectively. 



2. Group overlap for fixed probability quantiles 

In order to further study the differences in the performance 
of the considered measures, we apply another complimentary 
statistical test. Specifically, we take the a-quantile Q Q (Ai) of 
the distribution of Ai that corresponds to a given value A* (i.e., 
a ~ ri\jn) and consider a related decomposition of the (c, a) 
parameter plane based on the same quantile for the recurrence- 
based measures, i.e., 

S[(Q a (x)) := {(c, a)\x(c,a) < Q a {x)} 
S' 2 {Q a (x)) := {(c, a)\x(c,a) > Q a {x)} 

with x G {1 - DET, 1 - Lmean, 1 - C, 1 - £}. Since 
the quantile a has been kept fixed here, S[ and S' 2 contain 
the same numbers of elements (m and n 2 , respectively) as Si 
and S 2 - Hence, we are able to quantitatively assess the coin- 
cidence between the grouping based on Ai and x by consider- 
ing the relative frequency pf of (c, a) pairs that do not belong 
to the same group based on the two different measures. Fig- 
ure [TZt ) demonstrates that with respect to this criterion, the 
average path length C again shows (on average) the lowest 



An even more detailed characterization of classification er- 
rors is obtained in terms of the receiver operating character- 
istics (ROC) Ir33tl . In the ROC analysis, we compare the dis- 
crimination (Si, S'2) of the set of parameters (c, a) based on 
a fixed value of A* with another grouping (S[, S' 2 ) based on a 
variable threshold x* of the observable x, which now replaces 
the quantile Q a (x), The probabilities of correct as well as 

(v) (v) 

false detections of periodic behavior, pf' and pj ', respec- 
tively, are given as 

pW(A*,z*) = ISjnSil/ISil 

P f\x*,x*) - \s[ns 2 \/\s 2 \ 

where \S\ represents the cardinality (i.e., the number of ele- 
ments) of the set S. Varying x* over the full range of possi- 
ble values with A* simultaneously being kept fixed, we obtain 
a continuous curve in the {pf ,Pc )-plane, the ROC curve, 
which illustrates the trade-off between a high probability of 
correct and a low probability of false detections of periodic 
behavior (Fig. [T3K). The area under this curve, AU C, can 
be (among other statistics) used for quantitatively characteriz- 
ing the classification performance of different measures ll63ll 
(Fig. [T3b). Specifically, high values of AU C correspond to a 
low probability of false classifications (i.e., high specificity) 
and, simultaneously, a high probability of correct classifica- 
tions (i.e., high sensitivity). In contrast to all other kinds of 
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statistical tests, the ROC analysis suggests that in the specific 
setting studied in this work, among the four considered mea- 
sures, C is the most suitable statistics for discriminating be- 
tween periodic and chaotic behavior of the Rossler system, 
followed by C, DET, and Lmean- 
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FIG. 13: (A) ROC curves for A* = 0.01, and (B) area under the 
ROC curve (AUG) in dependence on A* for all four measures. 
For A* = 0.01, AUG takes the values 0.9279 (DET), 0.9090 
(Lmean), 0.9487 (C), and 0.9442 (£), respectively. 



One has to note that the overall classification error 

Pf = p f ( A* , Q a (.*) ) + pf ( A* , Q a (x) ) 
= p < f\\*,Q a (x)) + (l-pM(\*,Q a (x)) 



(13) 



(with pf (A* , Q a (x) ) = | S' 2 n Si | / 1 Si being the relative fre- 
quency of false detections of chaotic dynamics) still reaches 
values of around 10% even for the best performing measures 
(Fig. [l~2"t)). The main reason for this performance failure is 
the different sensitivity of Ai and the recurrence-based mea- 
sures close to the bifurcation lines (see Fig. H3V Apart from 
this, all four measures allow recovering the overall structures 
in parameter space comparably well as Ai. The obvious dif- 
ferences in the transition regions could be related to the use of 
short time series in regions with a complex bifurcation sce- 
nario, which are probably affected by problematic features 
such as longer transients before the attractor is reached (see 
Sec. Ill Bl i. different bifurcation scenarios characterizing the 
transition between periodic and chaotic dynamics across dif- 
ferent boundaries of the same shrimp, intermittency, or dif- 
ferent stiffness properties of the considered trajectories J3l0]. 



V. CONCLUSIONS 

We have proposed to use nonlinear recurrence-based char- 
acteristics of time series for exploring the parameter space of 
complex systems. This is of special importance when dealing 
with experimental data. Specifically, for distinguishing peri- 
odic and chaotic dynamics, in the recent literature estimates 
of the maximum Lyapunov exponent Ai from the correspond- 
ing ODEs have most often been used, which allow resolv- 
ing the borders of shrimps in a satisfactory way (see Fig. |3). 
This specific approach works well if the associated variational 
equations are provided explicitly. However, if these equa- 
tions are not available (as in the case of experimental data), 
the numerical estimation of Ai is typically much more chal- 
lenging, especially when dealing with short time series l32ll . 
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FIG. 14: Discrimination errors (black dots) for the quantile-based 
groupings Si, S 2 , S[, and S' 2 for A* = 0.01 (see Sect. |IVC2) for 
(A) DET (p f = 0.0923), (B) L ME an (0.1106), (C) C (0.0954), 
and (D) £ (0.0899). 



The recurrence-based measures used in this work have the ad- 
vantage that they can be properly estimated from rather short 
time series, so that they could prove advantageous also in sit- 
uations where the available amount of data is not sufficient 
for obtaining reliable estimates of Lyapunov exponents. Us- 
ing recurrence-based methods instead of Lyapunov exponents 
therefore has great potentials for the automatized discrimi- 
nation between different types of dynamics in many appli- 
cations. Specifically, the alternative discriminatory statistics 
introduced in this work could become very helpful not only 
in evaluating, but also already in designing both experimental 
and numerical studies, since the requirements concerning the 
necessary amount of data can be much easier matched. As a 
perspective, we therefore expect that systematic application of 
the proposed methods will open new fields of applications, as 
it has already been the case with the introduction of RQA 13411 . 
However, there are still situations where the consideration of 
Lyapunov exponents is clearly superior for detecting periodic 
windows in comparison to the recurrence-based approach, in 
particular when the governing equations are explicitly known 
or the available observations of the considered systems are 
very long. It will remain a subject of future studies to investi- 
gate in more detail under which specific conditions which of 
the two approaches is favorable. 

Another traditional idea to distinguish periodic behavior 
from chaotic dynamics is to use a properly defined surface 
in phase space for a Poincare section, since a periodic tra- 
jectory has only a finite number of intersection points with 
this surface, while a chaotic one renders an infinite number of 
crossings. Note, however, that constructing a proper Poincare 
section for a given set of ODEs l64ll is much easier than it 
is the case when working with time series, where one most 
often relies on some interpolation techniques to find cross- 
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ings on the section, which obviously introduces noise-like ef- 
fects 13211 . Furthermore, there are no universal criteria for 
choosing Poincare sections when scanning a large parameter 
space (e.g., Fig. [3j since the shape and the orientation of the 
attractors typically vary for different parameter values. 

In this work, we have been particularly interested in the 
problem of detecting specific periodic windows in parameter 
space, so-called shrimps, which are characterized by a rich bi- 
furcation scenario and self-similarity. Specifically, we have 
addressed the problem of numerically detecting shrimps in 
systems of ODEs, while the related question of the associated 
bifurcation scenarios remains a subject for future investiga- 
tions. We expect the recurrence-based methods proposed in 
this work to be particularly useful for this problem, especially 
concerning the properties of secondary shrimps and the quan- 
titative analysis of possible transient behavior. 

For properly identifying shrimps in parameter space based 
on recurrence plots, both measures from RQA and complex 
network theory have been applied. For the (c, a) parameter 
plane of the Rossler system, especially the recently proposed 
application of network measures to the recurrence matrix of 
complex systems l36l l39ll yields results coinciding well with 
those obtained using Ai. The used RQA measures in this pa- 
per need two parameters recurrence rate RR and minimum 
line length l m i n , RN measures depend exclusively on RR. 
We have to emphasize that this conceptual advantage comes 
at the cost of higher computational demands, especially when 
considering the average path length C of the RNs. How- 
ever, the clustering coefficient C, which has been found to 
perform best in our example when considering the most so- 
phisticated statistical evaluation (ROC analysis), can be cal- 
culated at significantly lower costs. Generally, C seems to 
be very well suited for distinguishing periodic from chaotic 
dynamics, even for time series sampled at a very high rate, 
such as the examples presented in this work. However, for 
high sampling rates, DET and Lmean can have a reduced 
discriminatory performance since the typical line structures 
become too similar 065I1 . We additionally note that both net- 
work measures have shown a slightly better discriminatory 
power for secondary shrimps than the two RQA measures 
DET and Lmean, at least for the specific setting used in 
this study. As we have not explicitly shown here, even larger 
performance errors can be observed for other line-based RQA 
measures. A more detailed investigation of possible impacts 
of other choices of l m i n and RR has not yet been systemat- 
ically performed, but will be subject of future studies. As a 
preliminary result, repeating all presented calculations with a 
different value of RR = 0.01 supports all results discussed in 
detail for RR = 0.02 in this paper. 

The two specific network-theoretic quantities considered 
here have already been shown to be able to detect dynami- 
cal transitions in both model systems and real-world time se- 
ries lHU. We note, however, that since both measures are in- 
variant under permutations of the time coordinate ll39ll . they 
exclusively capture the geometry of states associated to a spe- 
cific trajectory in phase space. This is a distinctive difference 
to most existing methods of time series analysis, which are 



related to the study of temporal correlations. We hypothesize 
that the good performance of network measures in detecting 
shrimps could be related to this fact, since network measures 
take only information about spatial correlations (i.e., neigh- 
borhood relationships in phase space) into account. 

While a rigorous theory interrelating the phase space prop- 
erties captured by RN measures to traditional dynamical in- 
variants is still missing, the algorithm for estimating the corre- 
lation entropy K2 from RPs (which has also been successfully 
applied for detecting shrimps in continuous systems lfl5ll ) has 
been justified theoretically rf35tl . Since the estimation of K2, 
however, involves several steps lfl5ll with some subjective is- 
sues (in particular, the choice of the scaling region for con- 
vergence), its computational demands are significantly higher 
than those of the complex network as well as RQA measures 
applied in this work. Particularly, a large range of values of 
RR (1% ~ 99%) is often considered to yield a well defined 
plateau of K2(RR), while the method proposed in this paper 
works with just one chosen value of RR, dramatically reduc- 
ing computational costs. Systematically calculating RN statis- 
tics as well as RQA measures can be easily automatized using 
a fixed recurrence rate RR, ideally with RR < 0.05 JHHHt]. 
Hence, we conclude that our approach allows a systematic 
numerical discrimination between periodic and chaotic dy- 
namics of a continuous system, which is more practicable 
than other possible techniques especially when systematically 
studying higher-dimensional parameter spaces. 

In relation with the problem of discriminating between two 
qualitatively different types of dynamics in a binary way, 
which has been discussed in this work, Lyapunov exponents 
have also found wide use in quantitatively characterizing the 
"chaoticity" of complex dynamics. Following the results from 
Tab. HU we emphasize that the recurrence-based characteris- 
tics considered in this work may be used for similar purposes. 
In contrast to the discriminatory power, our corresponding 
initial results suggest that the applied RQA measures may 
be somewhat better suited for this purpose than the network- 
based concepts. However, further detailed statistical analysis 
is necessary in order to provide further evidence that this re- 
sult holds in general. A detailed treatment of this question 
therefore remains subject of future studies. 
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